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ABSTRACT 

Analytical studies based on perturbative theory have shown that the moments of the 
Probability Distribution Function (PDF) of the local smoothed velocity divergence are 
expected to have a very specific dependence on the density parameter fi in the quasi- 
linear regime. This dependence is particularly interesting as it does not involve the 
possible bias between the galaxy spatial distribution and the underlying mass distri- 
bution. This implies a new and promising method for determining a bias-independent 
value of f2 based on a reliable determination of the velocity divergence PDF. 

In this paper we study the f2 dependence of the velocity divergence PDF and 
its first moments in a set of N-body simulations, using the so-called Voronoi and 
Delaunay methods. We show that this dependence is in agreement with the theoretical 
prediction, even while the number density of velocity field tracers has been diluted to 
a value comparable to that available in current galaxy catalogues. 

In addition, we demonstrate that a sufficiently reliable determination of these 
statistical quantities is also possible when the measurement of the galaxy peculiar 
velocities is restricted to the one component along the line-of-sight. Under ideal, noise- 
free circumstances we can successfully discriminate between low and high fi. 

Key words: Cosmology: theory - large-scale structure of the Universe - Methods: 
numerical - statistical 



1 INTRODUCTION 

The study of the cosmic velocity field is a very promising 
and crucial area for the understanding of large-scale struc- 
ture formation. Since the early work of Rubin et al. (1976) 
and Burstein et al. (1987) a lot of effort has been invested 
in the measurement of the large-scale velocity flows (see 
Dekel 1994 and Strauss & Willick 1995 for recent reviews 
of the subject). Particularly important for the analysis of 
these measured velocity fields has been the development 
of the parameter-free POTENT method by Bertschinger & 
Dekel (1989). Based on the plausible assumption of poten- 
tial flow it enabled the construction and study of the full 
three-dimensional velocity field in a fair fraction of the lo- 
cal Universe out of the measurements of galaxy line-of-sight 
peculiar velocities. This opened up and triggered a host of 
studies addressing various issues and aspects of cosmic ve- 
locity flows and provided a versatile ground for testing the 
scenario of large-scale structure formation. 

The cosmic velocity field is particularly interesting be- 
cause of its close and direct relation to the underlying field 
of mass fluctuations. Indeed, on these large scales the accel- 
eration, and therefore the velocity, of any object is expected 



to have an exclusively gravitational origin so that it should 
be independent of its nature, whether it concerns a dark 
matter particle or a bright galaxy. Moreover, the linear the- 
ory of the generic gravitational instability scenario predicts 
that at every location in the Universe the local velocity is 
related to the local acceleration, and hence the local mass 
density fluctuation field, through the same universal func- 
tion /(fi) « Q. os (Peebles 1980). As the linear theory pro- 
vides a good description on those large scales the use of this 
straightforward relation implies the possibility of a simple 
inversion of the measured velocity field into a field that is 
directly proportional to the field of local mass density fluc- 
tuations 8 = pi (p) — 1. Such a procedure can be used to 
infer the value of SI, through a comparison of the resulting 
field with the field of mass density fluctuations in the same 
region. However, the determination of this mass density fluc- 
tuation field through the measurement of the local galaxy 
density fluctuation field, 5 g , may be contrived. The galaxy 
distribution may be representing a biased view of the under- 
lying mass density fluctuation field. A common and rather 
simplistic assumption is that 8 g and S are related via a linear 
bias factor 6, 
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Fig ure 1. The PDF of the velocity divergence as given by Eq. 
(|10|) showing its dependence with fl and a. 



S g (r)=b5(r). 



(1) 



However, although several physical mechanisms have been 
invoked to explain such a linear bias model (see e.g. Dekel & 
Rees 1987), by lack of a complete and self-consistent theory 
of galaxy formation it should as yet only be considered as 
a numerical factor roughly describing the contrast of galaxy 
density fluctuations with respect to the mass density fluctu- 
ations. 

The comparison between the observed local galaxy den- 
sity fluctuation field and the local cosmic velocity field, in- 
voking equation ([[]), will therefore provide an estimate of 
the ratio, 



(2) 



Various studies, most notably the ones based on a compar- 
ison of the galaxy density field inferred from the IRAS red- 
shift survey of Strauss et al. (1990) and the local velocity 
field reconstructed by the POTENT algorithm (Bertschinger 
et al. 1990; Dekel, Bertschinger & Faber 1990; see Dekel 
1994), have yielded estimates of /3 in the range fj « 0.5— 1.2 
(see Dekel 1994, Strauss & Willick 1995, for compilations of 
results) . 

It is then of crucial interest to find ways to disentangle 
the contribution of SI and b to /3. A variety of methods using 
intrinsic properties of the large-scale velocity field have been 
proposed to achieve this. One such attempt is based on the 
reconstruction of the initial density field from the observed 
distribution of matter through the use of the Zel'dovich ap- 
proximation. The further assumption of Gaussianity of the 
initial density probability distribution leads to a constraint 
on Q (Nusser & Dekel 1993). Although promising, the quan- 
titative results of this approach may be questionable, as the 
Zel'dovich approximation provides non-exact results for the 
induced non-Gaussian properties of the density and veloc- 
ity field (Bouchet et al 1992, Bernardeau 1994a). Another 
interesting attempt is the one by Dekel & Rees (1994), who 
exploit the simple observation that voids have a maximal 



"emptiness" of S = — 1. On the basis of the corresponding 
analysis of the measured velocity field in and around a void 
region these authors inferred a lower limit on £1 of about 0.3. 

In this paper we focus on a method to determine the 
value of fl that finds its origin in a statistical analysis of 
the velocity field. The foundation for this method is formed 
by analytical work within the context of the perturbation 
theory for the evolution of density and velocity fluctuations. 
In particular, it focuses on the statistical properties of the 
divergence of the locally smoothed velocity field 8, which is 
defined as 



(3) 



H 



where v is the time derivative of the comoving coordinate 
x, and V = d/dx. The method, proposed by Bernardeau 
(1994a) and Bernardeau et al. (1995), exploits the relations 
between the lower order moments of the probability distri- 
bution function (PDF) of 6, and the explicit dependence on 
Q. of these relations. The specific form of these relations were 
derived under the assumption of Gaussian initial conditions. 

The viability of the method was demonstrated by 
Bernardeau et al. (1995), who used the strong ^-dependence 
of the skewness factor of 9, or rather the third normalized 
moment T3, 



T s = (e a )/(e 2 ) «n 



(J) 



to estimate successfully the density parameter in N-body 
simulations of structure formation. A tentative application 
of this method to the observed velocity field as processed by 
the POTENT method yielded results consistent with Q = 1. 
Moreover, subsequent work by Bernardeau & van de Wey- 
gaert (1996) showed that the theoretical predictions con- 
cerning the complete overall shape of the PDF of 6 were 
valid. They demonstrated this by comparison of the analyt- 
ically predicted PDF with the PDF determined from a large 
CDM N-body simulation with fl = 1. In order to deal with 
the major complication of obtaining clean and straightfor- 
ward numerical estimates of the statistical properties of the 
velocity divergence field from the discretely sampled velocity 
field, they developed two new techniques. These techniques 
exploit the minimum triangulation properties of Voronoi and 
Delaunay tessellations (Voronoi 1908, Delaunay 1934, see 
Van de Weygaert 1991, 1994 for references and applications). 

A Voronoi tessellation of a set of particles is a space 
filling network of polyhedral cells, with each cell being de- 
fined by one of the particles as its nucleus and delimiting the 
part of space closer to this nucleus than to any other of the 
particles. Closely related to the Voronoi tessellation is the 
Delaunay tessellation, a space filling lattice of tetrahedra (in 
three dimensions). Each of the tetrahedra in the Delaunay 
tessellation has four particles of the set as its vertices, such 
that the corresponding circumscribing sphere does not have 
any other particle inside. Through the duality relation of 
the Voronoi and the Delaunay tessellation it is possible to 
obtain one from the other. 

In the method based on the Voronoi tessellation - the 
"Voronoi method" - the velocity field is defined to be uni- 
form within each Voronoi polyhedron, with the velocity at 
every location within each cell being equal to that of its 
nucleus. The obvious implication of such an interpolation 
scheme is that the only regions of space where the velocity 
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divergence, as well as the shear and vorticity, acquires a non- 
zero value is in the polygonal walls that separate the cells. 
While the Voronoi method can be regarded as a zeroth order 
interpolation scheme, yielding a discontinuous velocity field, 
the "Delaunay method" can be seen as the corresponding 
first order scheme. Basically it constructs the velocity field 
within each Delaunay tetrahedron through linear interpo- 
lation between the velocities of the four defining particles. 
Evidently, the velocity field constructed by the Delaunay 
method is a field of uniform velocity field gradients within 
each Delaunay cell. In the following we consider the statis- 
tical properties of the local velocity divergence when it is 
filtered by a top-hat window function following a volume 
weighting prescription. In practice the local filtered diver- 
gence is computed on 50 3 grid points in all cases. Depending 
on the method used to define the velocity field, the filtered 
divergence is given by a sum of intersections of a sphere with 
either planar polygonal walls or tetrahedra (each of them be- 
ing multiplied by the local divergence) divided by the vol- 
ume of the sphere. In such a volume weighting scheme the 
filtered quantities do not depend on the local number den- 
sity of tracers. But of course the more numerous they are, 
the more accurately determined are the local divergences. 
For an extensive description and preliminary tests of the 
two techniques we refer to Bernardeau & van de Weygaert 
(1996). 

With the intention of demonstrating the potential and 
practical applicability of our method, this paper presents 
the results of a systematic study of the Q dependence of 
the moments and the PDF of the velocity divergence 8 in 
N-body simulations of structure formation, using the numer- 
ical schemes of the Voronoi and the Delaunay method. To 
this end, we will first recall the relevant theoretical results 
on the statistical properties of 8 in Section ^| thereby un- 
derlining the main features of importance to our study. This 
theoretical groundwork is followed by Section |§| containing 
the presentation of the numerical results of the statistical 
analysis of PM N-body simulations of structure formatio n in 
Q — 1 and Q < 1 universes. While the first subsection, § 3.1 
concerns the statistical quantities that were determined in 
a very large sample and thus with maximum attainable ac- 
curacy, the second subsection § 3.2 is devoted to the issue 



to what extent these results get affected when the sample 
is strongly diluted. The latter is particularly important as 
we are interested in the reliability of our method for sam- 
ples with a number density comparable to that of available 
galaxy catalogues. Also of immediate relevance for a practi- 
cal application is the question of how far the results of the 
statistical analysis get influenced when the velocity of the 
particles in the sample are known along only one direction. 
This issue, of crucial importance within the context of the 
statistical analysis of obs ervational catalogues, is treated in a 
third subsection, in § 3.3. Finally, following the successful ap- 



plication of our method under the circumstances described 
above, we conclude with a summary and a discussion of pos- 
sible complications and prospects for our statistical method 
to infer a bias-independent value of Q. 



2 PERTURBATION THEORY OF 

STRUCTURE FORMATION AND THE 
VELOCITY FIELD PROBABILITY 
DISTRIBUTION FUNCTION 

Perturbation Theory (PT) is extremely useful for the study 
and analytical description of the mildly non linear evolu- 
tion of density and velocity fields. In particular within the 
context of structure forming out of Gaussian initial condi- 
tions perturbation theory has been extensively developed 
in a large body of work (see e.g. Bernardeau 1994a, b). In 
the case of these Gaussian initial conditions the complete 
set of moments of the smoothed velocity and density fields 
can be computed analytically, in particular if these fields are 
top-hat filtered. The corresponding PDF can be computed 
through re-summation of the series of moments. 

One of the most straightforward and useful results in 
the context of perturbation theory is the relation between 
the third moment (# 3 ) and the second moment (8 2 ^ of the 
probability distribution function of 8, 



(e 3 ) = T 3 (e 2 ) 2 = n4. 



(•») 



The coefficient T3 depends on the cosmological parameter 
Q, on the shape of the power spectrum, on the geometry of 
the window function that has been used to filter the velocity 
field and even on the value of the cosmological constant A, 
although the latter is an almost negligible weak dependence. 
In fact, the dependence of T3 on SI is substantially stronger 
than the one of the equivalent coefficient for density field. 
For instance, for a top-hat window function and a power law 
initial power spectrum of index n, i.e. 

2\ 



P(k) = (S(k) 2 ) oc <6»(k) 2 ) oc k n , 
one obtains the following expression for T3 
-1 



T 3 



SI - 6 



f-(3 + n) 



(6) 



(7) 



As T3 can be directly determined from observations, one can 
use its strong dependence on Q to obtain an estimate of Q, 
as has been done by Bernardeau et al. (1995). 

More generally, Perturbation Theory enables one to in- 
fer the whole set of the cumulants (8 P ) to their leading order. 
All of them are related to the second moment via the rela- 
tion 



{en=T p (e 2 ) p -\ 



(8) 



and as in the case of T3 all the coefficients T v possess a 
strong dependence on the value of Q (Bernardeau 1994b). 
To a good approximation, this dependence on cosmological 
parameters can be written as 

r„(n,A) w Q(p } 2)0S T P (Q = 1, A = 0). (9) 

This property, given here for the moments, naturally ex- 
tends itself to the shape of the complete velocity divergence 
PDF p(8). This can be directly appreciated from the work 
by Bernardeau (1994b), who showed that the PDF can be 
calculated from its moments T p and the value of a 2 through 
a Laplace transform of its generating function ipg 



2%\o-', 



exp 



„2 "+" _2 



d8 . (10) 



The moment generating function ipg(p,,y), given by 
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-y) p 



(ii) 



p=2 



can be related to the spherical collapse dynamics in the cos- 
mology under consideration (see Appendix A). Although 
this calculation is almost intractable in the general case, 
it has been shown that at least in one particular case one 
can evaluate this expression analytically. In the specific case 
of a power law density perturbation spectrum with an in- 
dex n — — 1 in combination with a top-hat smoothing, it 
is possible to invoke some approximations that enable the 



derivation of a simple analytic fit for the PDF 
Appendix A) . This fit is given by 



(see 



p(0)d8 



([2k-1]/ K 1 / 2 + [A-1]/A 1 < /2 )- 3 '' 2 



exp 



2\a 



K 3/4(27r)l/2 ffe 

dO. 



(12) 



with k = 1 + 9 2 / (9\fl 12 ) , and A = 1 — 29/ (3fi ' 6 ) . 

The behaviour of this function p(9) has been illustrated 
in figure |l]for various values of Q, and og. Qualitatively, one 
can see that the dependence of the shape on fi reveals itself 
in two ways: (1) the location of its cut-off at high positive 
values of 9 and (2) the location of its peak. 

As for (1), the maximum value that 9 can obtain is 
known exactly and is not dependent on the approximations 
that have been invoked to derive expression ( |l2] ) 

max = 1.5 n & . (13) 

The value of # max determines the location of the cut-off, and 
therefore the maximum expansion rate in voids. The value 
of 1.5 is the difference in value of the Hubble parameter in 
an empty, il = 0, Universe and that of an Einstein-de Sitter 
Universe, fl = 1. Evidently, this is reflecting the fact that 
the interior of the deepest voids locally mimic the behaviour 
of an Q = Universe. Recall that the suggestion by Dekel 
& Rees (1994) of using the maximum emptiness of voids to 
constrain f2 is also based on a similar feature. 

Also quite sensitive to the value of fl are the position 
of the peak of the distribution function p(9), i.e. the most 
likely value of 9, and the overall shape of p(9). Using the 
Edgeworth expansion (Juszkiewicz et al. 1995, Bernardeau 
& Kofman 1995) one can show that the value of 9 for which 
the distribution reaches its maximum is given by 



^pcak 



T 3 1 



(14) 



In fact, a procedure exploiting this dependence of shape and 
peak location of p(9) will probably yield a more robust mea- 
sure of f2 than the maximum value of 6> max as it will be less 
bothered by the noise in the tails. 



3 THE Q DEPENDENCE IN NUMERICAL 
SIMULATIONS 

By means of numerical simulations we have investigated the 
discussed dependence of the PDF of 9. These N-body sim- 

* This fit is obtained through approximations which tend to lower 
the values of the moments (appendix A). The PDF ( jl2| ) presented 
here is actually more accurate if n m —1.3. 



ulations use a Particle-Mesh (PM) code (Moutarde et al. 
1991) with a 256 3 grid to follow the evolution of a system 
of 256 3 particles. For our project we used two simulations, 
one with Q having a value of Q = 1 and the second one 
of f2 < 1. By analyzing the latter at different time-steps 
we explore situations for different values of fi. The particle 
distribution in the two simulations corresponds to a density 
and velocity fluctuation field with a P(k) oc fe _1 spectrum. 

As can be seen in Table ^ the variances ae do not differ 
significantly for the different values of Q for a given filtering 
radius. The fact that the values of the variance are com- 
parable simplifies a comparison of the PDF substantially, 
which makes the interpretation in terms of the intrinsic f2 
dependence more straightforward. 

3.1 Measurements with a large number of tracers 

The first step of our analysis concerns an exploration of the 
velocity field using a large number of tracers. For this study 
the number of selected particles in each simulation is about 
70,000, which for a cell radius of about 6% of the box size 
leads to a mean number of 67 particles per cell. The selec- 
tion procedure used here is deliberately biased towards low- 
density regions by inducing it to retain a uniform density 
of particles all over the simulation box. Except for its goal 
of achieving a better velocity field coverage of low-density 
regions such a selection bias is not expected to influence the 
velocity field analysis. 

The methods that we use to analyze the simulations are 
exactly the same ones as described by Bernardeau & van de 
Weygaert (1996). In fact, at this stage we only used the 
Delaunay method to calculate numerically the shape of the 
PDF of the velocity divergence. The results for the case of 
a large number of velocity field tracers are shown in Fig. ^| 
The results are in good agreement with the theoretical pre- 
dictions for the values of T3 and T4 (see Table 0), as well as 
with the theoretical shape of the PDF (Fig. [l]) . 

It is in particular worth noting that the specific features 
expected from equation (|l2] ) are indeed confirmed by the nu- 
merical results. Notably, the locations of the cut-off, which 
are very sensitive to rare event discrepancies, are well repro- 
duced (solid and dotted lines) . Moreover, as can be observed 
from the insets, also the position and shape of the peak have 
been reproduced very well (solid lines), providing a strong 
discriminatory tool between different values of Q. 

Within the context of these observations, we should is- 
sue a few side remarks. Although the shape (h_2|) is very at- 
tractive because it is a close analytic form, one should have 
in mind that it is only approximate. Indeed it is derived 
from an approximate expression for the cumulant generating 
function. The differences do not reveal for the overall shape 
(the logarithmic plots) but are significant for the shape of 
the peaks. For calculating the theoretical predictions we are 
then forced to use a more accurate description of the cu- 
mulants. To achieve this we use the relations (A. 13, A. 14) 
with n = —0.7 (the expression [l^ corresponds to n = — 1) in 
the integral (A. 15) which is then computed numerically. It 
is still an approximate expression, but it yields the correct 
value for T3 and a very good approximation for the higher 
order cumulants. We should emphasize that this slight mod- 
ification is only instrumental in obtaining the correct shape 
of the PDF around its maximum, which is indeed almost 
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Figure 2. The PDF of the velocity divergence for various values of Q. The dotted lines correspond to the approximate analytic fit 
(Eq.p^) and the solid lines to the theoretical predictions using (A. 15) and (A. 14) with n = —0.7 obtained for the measured values of a 
and Q. The dashed lines are the predictions for Q = 1 and the same variance. The numerical estimations have been obtained using the 
Delaunay method. 



entirely determined by the values of the low order moments. 
This may be understood for example from the properties of 
the Edgeworth expansion, for which we refer to Juszkiewicz 
et al. 1995 and Bernardeau & Kofman 1995. 

3.2 The effects of dilution 

In order to check the robustness of the results when only a 
limited number of tracers for the velocity field is available, 
we performed numerical experiments where only 10,000 par- 
ticles are used to trace the velocity field. The selection of the 
sample points in this diluted sample is completely random 
and does not invoke the specific biased selection procedure 



that was used in the case described in the former subsection. 
For this case of diluted samples, we used both the Delaunay 
and the Voronoi methods for analysis. 

Figure ^ shows the PDFs obtained with both meth- 
ods, for an £1 = 1 simulation and for an fl = 0.33 simula- 
tion. Both the Voronoi and the Delaunay methods appear 
to yield numerical results that are in reasonably good agree- 
ment with the predicted PDF. Particularly encouraging is 
the result born out by the insets, namely the fact that the 
shape of the peak can still be used as a strong discriminatory 
tool between different values of £1. 

Investigating Fig. ^ in somewhat more detail, we can 
observe that the results obtained by the two methods are af- 
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Figure 3. Effects of dilution for Q = 0.33 (left panels) and £7 = 1 (right panels). The PDF-s have been obtained with only 10,000 tracers 
in total which gives an average of about 10 particles per cell. The squares show the results of the Voronoi method and the triangles the 
results of the Delaunay method. The solid lines correspond to the theoretical predictions for the variance obtained with each method 
and the "right" assumption for Q and the dashed line with the "wrong" assumption for each case and method. 



fected somewhat differently by the dilution procedure. The 
PDFs obtained with the Voronoi method generally possess 
a less sharply defined tail at the side of the high 9 values. 
This effect appears to be stronger for the low Q. case. This 
behaviour may originate in the fact that the divergence is 
localized to a limited part of space, non-zero values being 
confined to the walls of the tessellation. In such a situation 
Poisson like errors in the measurements are expected to be- 
come particularly prominent in the heavily diluted areas of 
the void regions. An additional difference is a slight under- 
estimation of the values of the T p coefficients by the Voronoi 
method (see Table [j]). Unfortunately, we do not see a pos- 
sibility to correct for such effects. On the other hand, the 
Delaunay method seems to be more robust against such ef- 
fects. However, it tends to underestimate the value of the 
variance and the higher order moments. The latter is prob- 
ably a consequence of the fact that the effective filtering 
radius tends to be larger. 

3.3 The effects of reducing the information to one 
velocity component 

A major complication in the analysis of velocity fields un- 
der practical circumstances is the fact that only the velocity 
component along the line of sight can be measured. This 
therefore forms the second issue that we address in this pa- 
per. 

As yet we restrict ourselves to an artificial situation with 
ideal measurements, not yet to investigate the determination 
of the statistical quantities associated with the velocity field 
under realistic circumstances. To simplify furthermore our 
investigations we assume that we can use the approximation 
of an infinitely remote observer so that the radial velocity 
can be identified with one velocity component, namely the 
a>direction in the following. More specifically, we address the 
effects of the reduction of information concerning the veloc- 



ity field traced by a diluted sample of points. In principle, 
the fact that the velocity field is only known along one di- 
rection should not pose any problem. In the usual structure 
formation scenarios based on gravitational instability the 
large scale velocity field is expected to be non-rotational, 
implying it to be a potential flow and therefore the gradient 
of a potential that can be inferred from the measurement 
of only one component of the velocity (Bertschinger et al. 
1990, Dekel et al. 1990). 

In the Voronoi method non-zero values of the divergence 
8 are restricted to the walls of the Voronoi tessellation, where 
the local divergence is given by 



6> wa ii =n.Av, 



(15) 



with n being the normed vector orthogonal to the wall and 
Av the difference of the velocities on the opposite sides of 
the Voronoi wall. The expression for the local vorticity is 
given by 



<^wall 



n x Av. 



(16) 



Assuming potential flow, and hence a zero value of ui, this 
implies relations between the various components of Av and 
the following expressions for 6? wa n: 



Av x 



Av v 



Av z 



(17) 



In practice, this introduces the numerically unstable opera- 
tion of dividing by one component n as it can be arbitrarily 
close to zero. It may therefore be more reasonable to try 
to estimate the value of (9 wa n using the stable, but ad-hoc, 
prescription of 



gcstim. 
'wall 



Av x 



(18) 



where e is a small parameter of the order of e ~ 0.1. 

Note that such a prescription is not self-consistent in 
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6 8 



Figure 4. Effects of reduction of the information on the velocity field to only one component for Q = 0.33 (left panel) and f2 = 1 (right 
panel). We used the estimator ( |l8[ ) with e = 0.1 to determine numerically the local divergences. The solid lines correspond to the model 
Q obtained from (Q) with the parameters /i and <r e given in Table la The dashed lines are the resulting shapes of the PDF-s when 
one uses the same parameters but the "wrong" value for f2 (1 instead of 0.33 in the left panel, 0.33 instead of 1 in the right one). 




div Del , 67 pts/cell div Vor , 10 pts/cell, ID eps = 0.2 

Figure 5. Scatter plots that show the differences of the various estimators of the local divergence measured at 500 different random 
locations. 
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reproducing the full 3D velocity field .[jjndeed, depending on 
the way one goes from cell A to cell B - and there are in- 
finitely many ways to go from A to B, even while they are 
direct neighbours - one will not necessarily find the same 
value for #^ S aii m ' using the equivalent prescriptions for Av y 
or Av z . The method that we adopt here will therefore cer- 
tainly not be the most accurate method. However, we may 
expect it to be a reasonable approximation, and will there- 
fore use it here to illustrate the properties in which we are 
interested. 

Figure ^] displays the results obtained with ( |l8| ) for 
e = 0.1. Evidently, the results get affected to quite some 
extent by the transition from ( lfi|) to (|l8|). This leads to the 
key question whether it is sti 1 feasible to reliably recover 
the statistical information on 9 or not. As can be observed 
from the scatter plots in Fig. [|, the change from a situation 
in which one has knowledge of the full velocity to one where 
this has been limited to only one component thereof intro- 
duces a large scatter. However, we found that it is possible 
to define a meaningful representation of the scatter plots in 
terms of the following empirical description: 



= H (9 + e) . 



(19) 



Within this expression the coefficient jj, is a constant with 
some fixed value. The scatter is represented by the quantity 
e, a Gaussian random variable whose value is independent 
of 9 and which has a vanishing mean. Using the measured 
values of the variance and the skewness of the distribution 
it is possible to estimate the value of fi and the value of the 
rms fluctuation a e of e. From ( pj| ) one can readily infer that 



Cestii 



= M 



T 3 



4 



T 3 



(20) 



where ag is the exact rms fluctuations of 9 and T3 its third 
cumulant, while cr cst im. and T3 cs tim. are the corresponding 
estimated values. By solving this set of equations one can 
find the values of fi and a e . Their values for various Q and 
e are listed in Table [| 

Moreover, significant within the context of the ultimate 
goal of developing an unbiased estimator of Q, is that these 
parameters were found to be almost independent of the value 
of Q. More specifically, it turns out that the value of fj, only 
depends on the adopted value of e, whereas a e is indepen- 
dent of fl and only marginally dependent on e. This can 
clearly be appreciated from the bottom right-hand panel of 
Fig. n, which demonstrates that the two estimations of 9, 
one based on e = 0.1 and the other on e — 0.2, are basically 
proportional to each other. It may therefore be argued that 
it is quite natural to expect that the noise e introduced by 
this method is somehow intrinsic to the distribution. 

In Appendix B we describe an extremely simple model 
based on the assumption that the relative velocity of two 



T This can be readily appreciated from the fact that the normal 
n to a wall in the Voronoi tessellation is proportional to the vec- 
tor Ar between the two points on each side of the wall. Thus, 
if it were possible to build a consistent velocity filed from the 
constraints (17) it would imply Av oc Ar yielding not only a 
vanishing vorticity but also a vanishing shear. 



Table 2. Numerical values of the parameters /1 (bias) and a e 
(noise) introduced by the calculation (^) of the velocity diver- 



gence. 



tt = 0.33 



0.1 



0.2 



0.73 
0.25 



0.56 
0.27 



= 1 



0.1 



0.2 



0.75 
0.24 



0.57 
0.26 



Table 3. Analytical estimations of the parameters /i and cr e (see 
Appendix B). 



Q = 0.33 e = 0.1 e = 0.2 

f*(e) 0.60 0.49 

CTe(e,ag) 0.20 0.23 



U = 1 



e = 0.1 e = 0.2 



u{e) 0.60 0.49 

<r e (e,cr s ) 0.19 0.22 



particles is proportional to their relative position. This al- 
lows us to compute analytically the parameters /j, and a e 
entailed by the use of the numerical scheme (^). One can 
show that they are both independent on Q, and only depend 
on e and a$, with a e < ag. These analytical predictions are 
listed in Table []. They appear to be fairly close to their 
numerical measurements. The discrepancy between the an- 
alytical and numerical estimations of these parameters is due 
to the loss of information associated with the projection of 
the velocity from three to one dimensions. Although this is 
not fairly represented by our model, the discrepancy seems 
to be quite small. 

On the basis of the simple model described in Eq. ( |l9| ) 
it is possible to reconstruct the corresponding shape of the 
distribution of # es tim, 

,a \ f + °° Astim. s exp(-e 2 /2/cr^) de , 

Pcstim. (flcstim. = / P 6 , ^75 ~ 21 

where p{9) is given by Eq. (^). A comparison of the result- 
ing distribution ( |2l| ) with the measured histograms is shown 
in Fig. ^| The agreement appears to be quite good, rendering 
this phenomenological description a quite valuable one. 

Evidently, we are still able to clearly distinguish be- 
tween the scenario with a high value of f2 and the ones with 
a low value, although the distinction is not as clear as in 
the previous cases based on ideal sampling circumstances. 
Looking into some detail, we come to the conclusion that 
the signal has been diluted somehow by the noise e, and 
that there is a competition between the true rms of the di- 
vergence ag and the value of <r e . In this respect it is also 
important to note that it is crucial for a successful deter- 
mination of Q to have good estimates of both /i and a e . 
Once these parameters have been determined, and as long 
as ag > <7 e , we can see no further problem in distinguishing 
between the different cosmological scenarios. Fig. ^ gives an 
idea of the magnitude of the discrepancy that one gets when 
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Table 1. Cumulants from the Perturbation Theory and as estimated by the various numerical methods. 



^ tracers per cell 


cumulants 




J i — r\ nn 

1 £ — u.oo 


o — n a 


O — 1 


CO 

Perturbation Theory 


T 3 PT 
Tf T 


-4.5 
31.41 


-3.33 
17.22 


-2.97 
13.67 


-1.71 
4.55 


« 65 


CT Dcl 


0.37 ±0.01 


0.41 ±0.01 


0.42 ±0.01 


0.38 ±0.005 


Delaunay Method 




-4.52 ± 0.2 


-3.3 ±0.2 


-2.9 ± 0.2 


-1.54 ±0.1 




j>DeI 


31.5 ±2.5 


17.3 ±3 


12.8 ±2.5 


3.1 ±0.7 


» 10 


CT Dcl 




0.36 ±0.01 


— 


0.35 ±0.005 


Delaunay Method 


^3 




-3.8 ±0.2 




-1.84 ±0.1 




j>Del 




24.7 ± 3 




5.4 ±1 


a 10 


O-Vor 




0.40 ±0.01 




0.38 ±0.003 


Voronoi Method 


^Vor 




-3.2 ±0.2 




-1.51 ± 0.04 




T^Vor 




19.4 ±5.5 




3.62 ± 1.5 


» 10 
one velocity 
component, e = 0.1 


^Vor 
rpVot 

J 4 




0.35 ±0.006 
-2.40 ±0.3 
11.9 ±4 




0.34 ±0.009 
-1.14 ±0.07 
2.9 ±2.3 



one assumes a wrong value of Q. In the left panel this con- 
cerns the case wherein one take a value of Q = 1 instead of 
the actual value of Q = 0.33, while in the right-hand panel 
it concerns the reverse situation. 

At this point it is worthwhile to stress once more that so 
far we only tested the method for ideal measurements of the 
line-of-sight velocities. Noise in the measured values of these 
velocities was not yet taken into account. Moreover, besides 
the fact that this noise has quite a large value, an evalu- 
ation of its influence is substantially complicated as it not 
only concerns random measurement errors but also contains 
contributions from a plethora of, partially un-understood, 
systematic effects. In a work in preparation we will attempt 
to develop specific techniques for estimating the velocity di- 
vergence PDF from such noisy line-of-sight measurements. 



4 DISCUSSION AND CONCLUSIONS 

We have tested and confirmed the validity of the strong £7 
dependence of the Probability Density Function (PDF) of 
the velocity divergence that had been predicted on the basis 
of analytical Perturbation Theory calculations. These tests 
are based on a numerical analysis of N-body simulations of 
structure formation in a Universe with a Gaussian initial 
density and velocity fluctuation field. On the basis of this 
verification we may conclude that the analytical predictions 
of Perturbation Theory yield very accurate results for a wide 
range of cosmological models. 

The main practical implication of our work is the basis 
it offers for a potentially very valuable and promising estima- 
tor of the value of Q, an estimator independent of a possible 
bias between the distribution of galaxies and the underlying 
matter distribution. The successful tests presented in this 
work demonstrate the validity of the equations of Perturba- 
tion Theory that form a basis for the estimates of Q which 
are based on the statistical properties of the velocity diver- 
gence field. Their unbiased nature finds its origin in the fact 
that the relations between the various statistical moments 



do not contain any explicit dependence on a bias between 
the galaxy and the matter distribution. Moreover, the esti- 
mated values of f2 are even more direct and straightforward 
to interpret as the relevant statistical relations only involve 
a very weak dependence on the cosmological constant A is 
expected to be very weak (Bernardeau 1994a, b). 

Not only relations between the statistical moments of 
the 6 distribution, also the shape and general functional be- 
haviour provide a useful indicator for the value of Q. When 
we focus on the details of this functional behaviour of the ve- 
locity divergence distribution function - illustrated in Fig. ^ 
and |3|) - we can draw a few conclusions with regards the 
practical feasibility of obtaining reliable estimates from the 
shape of p(6). Both location and shape of the peak of this 
distribution appear to be robust indicators of the value of 
Q. On the other hand, the location of the maximum of the 
divergence 6 - i.e. the cutoff value of p(8) - appears to be 
much more sensitive to a poor sampling. This may make it 
harder for it to provide reliable estimates of fi from currently 
available observational catalogues (see Fig. ^). However, it 
is all the more encouraging that even on the basis of the 
cutoff value we obtained a reasonable agreement between 
theoretical predictions and numerical measurements. 

Finally, we addressed one further crucial issue towards 
an application of our estimation procedures to real data sets. 
This concerns the problem of not being able to obtain di- 
rectly the full three-dimensional velocity field. Instead, the 
velocity of a galaxy can only be measured along the line-of- 
sight. In a preliminary attempt to study the consequences 
of this fact for the feasibility of our method, we introduced 
a partially empirically denned extension of our method. De- 
spite of the extreme crudeness and rather ad-hoc nature of 
this algorithm to reconstruct the full velocity field, it is quite 
encouraging that we are able to distinguish between the ve- 
locity PDF obtained in a flat Universe and that obtained 
in an open Universe. The major obstacle towards a success- 
ful application of our methods therefore appears to be the 
one of noisy data sets and systematic sampling errors. We 
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have not yet dealt with these problems, deferring them to a 
forthcoming paper. 

An additional and useful application of our numerical 
work involves a test for structure indeed having emerged 
through the process of gravitational growth of an initially 
Gaussian random density and velocity field. Having shown 
Perturbation Theory to be valid, we can exploit its predic- 
tion that the PDF of 8 is only dependent on a few parame- 
ters, in particular ag and SI. If no values of ag and SI can be 
found to produce an acceptable fit to the observed velocity 
field, this will force us to conclude that it is unlikely that 
the structure developed as described within the standard 
framework of gravitational instability and Gaussian initial 
conditions. In this context it is interesting to point out that 
a negative skewness has been observed in the currently avail- 
able datasets (Bernardeau et al. 1995), which is an indication 
in favour of standard scenarios. 

Summarizing, we may conclude that the combined ma- 
chinery of the analytical perturbation theory results and the 
developed numerical methods and their application on the 
intrinsic statistical properties of the velocity field provides 
us with a reliable new estimator of the cosmological density 
parameter SI. This estimator is all the more useful as it is 
one of the very few which will yield values of SI completely 
independent of galaxy-density field biases and almost inde- 
pendent of the value of A. 
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APPENDIX A: CALCULATION OF THE 
PROBABILITY DISTRIBUTION FUNCTION 
OF 8 

In the case of a top-hat window function it is possible to 
evaluate the whole series of cumulants of 8 at their lead- 
ing order. This makes it feasible to construct, at least in 
principle, the complete PDF of the local top-hat smoothed 
velocity divergence field. However, in order to construct a 
convenient closed form of the PDF we to introduce some 
approximations. Here we will go through some of the ana- 
lytical results and applied approximations that were used in 
order to obtain simplified final expressions. 

A key element in the construction of the PDF of 8 is 
the moment generating function ipg and the close relation 
that was discovered to exist (see e.g. Bernardeau 1992) be- 
tween this function and the dynamics of spherical collapse 
in the background Friedmann-Robertson- Walker universe. 
More specifically, the generating function ipg(Sl,y), defined 
as 

<p e (Sl,y)=J2 -T P (fi) LJ f L , (Al) 

p=2 

where T p are the reduced cumulants of 8 introduced in (|J), 
is determined from a function Qg(Sl,r) whose behavior can 
be deduced from spherical collapse dynamics in the cosmol- 
ogy under consideration. These two functions ipg(Sl,y) and 
Qg(Sl,r) are related through the system of equations 

tpg{Sl,y) = yGe(Sl,T)-~yT-£-g e ((l,T), (A2) 

2 clr 

t = -V^e{Sl,T). (A3) 

For the solution of this system it is therefore necessary 
to first determine the explicit relationship of Qg{Sl,r) to 
the dynamics of spherical collapse. In order to accomplish 
this, it is useful to introduce the functions G$ c {r) and 
gl c [Sl(a),f(Sl, A)r]. The first one of these, Gf C {r), is de- 
fined to be the nonlinear density contrast of a spherical 
perturbation of initial over-density — r. For exact analyti- 
cal expressions for Qf c , which exist only if A = 0, we refer 
to appropriate textbooks. Directly related to the expression 
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for the density contrast is the function describing the local 
departure from the Hubble expansion, 



gi c [n(a),f(n,A)r] 



d n sc i \ . 
. a _g s (ajT)+ 



/(fl, A) ^f>, r)/ [1 + Gl C (a, r)] 



(A4) 



After having introduced the above two functions and 
after having written down their explicit expressions, the 
generating function Qb{t) in equation ( p"l| ) can be obtained 
through the solution of the system of equations (Bernardeau 
1994b), 



&(/(«. A)t) 
&(t) 



/(n,A)r 



^[H-C^t)] 1 "/*,) 



<T([l + g^(r)] 1/3 ^o) 
<r(#o) 



(A5) 



In principle, the above system of equations can be fully 
solved to get full-fledged expressions for the coefficients T p , 
like e.g. the ones for T3 and T4 



26 

Ts(n = l,A=0) = -y-7i, 



(A6) 



T 4 (n = l,A = 0) = 



12088 338 



441 



where the 7 P are the successive logarithmic derivatives of 
the variance with the smoothing scale, 



1p 



dHogqg(fl) 
dlog p i? 



(A8) 



For practical purposes, however, it is preferable to define a 
simplified set of equations that forms a reasonable approxi- 
mation to the original system. Several useful approximations 
can be applied. The first one is to assume that the power 
spectrum of the density and velocity fluctuations can be ac- 
curately described by a power law, 



P(k)<xk n , a 2 (R) oc _R~ (n+3) , 



(A9) 



leading to 71 = — (n + 3) and 72 = in Eqs. (A6) and (A7). 
Although this may not be exact in general, it may be justi- 
fied by the fact that we saw previously that the corrections 
in the values of T3 and T4 induced by a variation of the index 
n enter only weakly. 

A second approximation exploits the fact that the fl 
dependence of Qf c (t) is extremely weak (see Bouche t et 
al. 1992, Bernardeau 1992). This implies that equation (A4) 
reduces to 



gl c [f(n,A)r] 



-n°- 6 T Agf °( T )/ [1 + gf c ( r )] . (A10) 



Moreover, it is therefore also reasonable to approximate the 
function gf c (r), independent of the value of f2 involved, by 
the simple expression for gf c (t) for the case £2 = 0, 



arw«(i + £) 



2r\- 3 / 2 



- 1 . 



(All) 



so that equation (A4), via equation ( A-10 ), lead to the ap- 
proximate relationship 



gl c {T)~r(i + . 2r 



3S1 



(A12) 



Using the above app rox imate expressions in solving the sys- 
tem of equations in ( |A3[ ) then yields an approximate expres- 
sion for the relation between the spherical collapse density 
contrast r and the functions g$ c (t) and Q/j c {r), 

-) ( - +3)/6 [(i+^r 2/3 -i], (A13) 



a + g s 



and 



r 1- 



2g e \ ( " +3)/6 

3/(fi) J 



1 - 



:»)) 



2g B 



(A14) 



While reasonable, these approximations turn out not to be 
as good for the velocity field as for the d ensity field . For 
example, using the approximations (All) and (A12) one 
would obtain S3 = ^<5 3 ^ / (^S 2 ^ — 5 — (n + 3) instead of 
S3 = 34/7 - (n + 3), and T 3 = -4 + (n + 3) instead of 
T3 — —26/7 + (n + 3). This means that for example for 
n = — 1 the relative errors in S3 and T3 are in the order of 
5% and 15%. 

While the above expose explains the strategy for cal- 
culating the coefficients T p to any order, the major purpose 
of calculating the generating function (|ll|) is to construct 
the PDF of the local velocity divergence. This is achieved 
through a Laplace inverse transform (Balian & Schaeffer 
1989), 



p(0)d8 - 



exp 



<fie(y) 



+ 



y6 
„2 



(A15) 



where a$ is the variance of the distribution. The saddle point 
approximation of this integral 

1/2 

I I — -r( , i -r \ / ( y I ~ \ 

Pe(e)d6 



1 - rg e (r)/g' e (r) 



2-koI 



xexp ( ) d9, g e {r) =1 



(A16) 



yields an accurate prescription for the shape of PDF, in 
particular around its maximum (see Bernardeau 1994b for 
more details). 

B y su bseque ntly invoking the approximate expres- 
sions (|All[) and ([A12j) , and using the fact that then equa- 
tions QA13j ) and ( |A14[) can be inverted, one can find simple 
analytic fits [Eq. (11 2h] for the special case of n = — 1. 



APPENDIX B: 6 FROM ONE VELOCITY 
COMPONENT: AN ILLUSTRATIVE EXAMPLE 

It is possible to estimate the bias fi and the variance of the 
noise cr e involved in our calculation of 9 [Eq. (g ^)] with 
the simple following model : let us assume that the difference 
of velocity between two given neighboring points is collinear 
to their relative position Ar 

Av = 6^, (Bl) 

leading to # wa ii = ©• This assumption is generally not true, 
even for an irrotational flow, because it overlooks any pos- 
sible shear. However this should be enough for a rough esti- 
mation of n and cr e , and to show that they do not depend di- 
rectly on Q. Considering that only the velocity along the di- 
rection x is known, our estimate of the divergence [Eq. (|l8|)] 
yields 
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= Av x = -£2- . (B2) 

ni + e raj + e 

If we now assume that the direction of Ar is randomly dis- 
tributed in a 3 dimensional space, the bias is, according to 
our empirical parametrisation (|l^), (the subscript "wall" of 
9 will be dropped hereafter) 

= (flctim. g) = / _n§_\ = j _ ^ arctan ( £ -im (B3) 
0^ \nz + e/ x ' 

This implies that the noise 

e = 6 ( . . - 1 ] (B4) 

has a vanishing average, and its variance is 
! / „4 



M \ (nl + e)* 



2 



,2 



1+ 2(TT^~¥ arctan ( e " 1/2 ) 



2 

/' 



05) 



. < 0.562 8 2 for < e < 1. 
We can now assume that is itself a random variable 
independent of n x such that (0) = (0) = and its vari- 
ance (0 2 ) = o" 2 is given in Table |l| This variance does not 
have an explicit dependence on fl, and therefore a e does 
not, while it does depend on the evolutionary stage of the 
system and the scale that is considered. On the other hand, 
H is clearly completely independent of Q. Analytical results 
for the bias and the noise are given in Table ^ and can be 
compared to their respective numerical values listed in Ta- 
ble § 

Taking this model at face value would tell us to put 
e = giving /i = 1 and e = 1, i.e. # es tim. would be a per- 
fect estimator of 9. This is of course not true because of the 
neglected shear component and of the vorticity likely to be 
present at small scale, not to mention the instability of such 
a numerical estimate. However, in spite of its extreme crude- 
ness, this model gives fairly good analytical predictions of 
the 'intrinsic' bias and noise introduced by our numerical 
estimates of the velocity divergence, and indicates that the 
remaining noise due to shear and vorticity is almost negligi- 
ble. 
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